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ABSTRACT 

Indirect detection of high-energy particles from dark matter interactions is a promis- 
ing avenue for learning more about dark matter, but is hampered by the frequent coin- 
cidence of high-energy astrophysical sources of such particles with putative high-density 
regions of dark matter. We calculate the boost factor and gamma-ray flux from dark 
matter associated with two shell-like caustics of luminous tidal debris recently discov- 
ered around the Andromeda galaxy, under the assumption that dark matter is its own 
super symmetric antiparticle. These shell features could be a good candidate for indirect 
detection of dark matter via gamma rays because they are located far from the primary 
confusion sources at the galaxy's center, and because the shapes of the shells indicate 
that most of the mass has piled up near apocenter. Using a numerical estimator specifi- 
cally calibrated to estimate densities in N-body representations with sharp features and 
a previously determined N-body model of the shells, we find that the largest boost fac- 
tors do occur in the shells but are only a few percent. We also find that the gamma-ray 
flux is an order of magnitude too low to be detected with Fermi for likely dark matter 
parameters, and about 2 orders of magnitude less than the signal that would have come 
from the dwarf galaxy that produces the shells in the N-body model. We further show 
that the radial density profiles and relative radial spacing of the shells, in either dark 
or luminous matter, is relatively insensitive to the details of the potential of the host 
galaxy but depends in a predictable way on the velocity dispersion of the progenitor 
galaxy. 

Subject headings: dark matter, galaxies: individual (M31), galaxies: kinematics and 
dynamics, gamma rays: galaxies, methods: analytical, methods: numerical 



1. Introduction 



The nature of the dark matter is one of the foremost questions in astrophysics. A lthough 
most astrophysicists agree that it is probably some kind of particle (jBertone et al.l l2005l ) , there 
are to date no conclusive detections. Countless experiments are attempting to do so, whether by 
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detecting dark matter particles directly ( Asztalos et al.ll2004l : lKinion et al.ll2005l : lGeralis et al.ll2009l: 



Ahmed et al.ll2009al lbl: lFenell20ld ). creating them in the laboratory (jBaltz et al 



Hooper k. Balta 120081) . or observing the st andard-model byproduct s of interactions b etween 



them (jChang et aljboosl : lAdriani et"alll2009al lbl: lAbdo et al.ll2009al . [2OI0I : IScott et alJboiol l. This 



last method, usually referred to as "indirect detection", usually assumes that the dark matter 
particle is its own antiparticle, so that an interacting pair of dark matter particles self-annihilates 
to produce various kinds of standard model particles. This assumption is motivated by predictions 
of supersymmetry that the dark matter could be the lightest supersymmetric partner (LSP) of a 
standard model particle, and by the cosmological result that such a particle with a mass of between 
20 and 500 GeV would have been produced in the early universe in sufficient number to resolve the 
discrepancy between the energy density of luminous matter observed today and the total energ y 
density of matter required to explain the gravitational history of the universe (jBertone et al.ll2005l ). 



Indirect detection, because it relies on observing the products of pairwise annihilation, has a 
signal strength that varies as the dark matter density squared. It is therefore most effective in re- 
gions with the very highest number density of dark matter particles. Because dark matter interacts 
with luminous matter primarily through gravity, these are often the same regions where the density 
of luminous matter is highest, such as the centers of galaxies (jBergstrom et al.lll998l ). Although 
the signal from pairwise dark matter annihilation may be highest in those regions, it also suffers 
from confusion with high-energy astrophysical sources like pulsars a nd X-ray binaries, which tend 
to be concentrated wherever the density of luminous matter is high (jAbdo et al.ll2009bl ). However, 
in some instances the dark matter density can be high while the luminous matter den sity is low; for 



exan iple, in the recently discovered ultra-faint dwarf galaxies orbiting the Milky Way (jStrigari et al 



20071 ). Since the dark matter has a lower kinematic temperature than luminous matter, it may also 
have more small-scale structure than luminous matter, increasing or "boosting" the production of 
standar d model particles by pairwise annihilation above the level predicted for a smooth distri- 
bution (Peirani et al. 2004: Diemand et al. 2007: Kuhlen et al. 2007: Mohavaee. Shandarin. Silk 



20071 : iBringmann et al.ll2008l : iPieri et al.ll2008l : IVogelsberger et al.ll2008l ). Cases where there is no 
confusion between particles produced by pairwise dark matter annihilation and those produced by 
high-energy astrophysical sources offer a high potential for a confirmed indirect detection, provided 
that the signal is still detectable. 

One possible scenario for indirect detection is the debris created by a merger between a 



larger host galaxy and a smaller progenitor galaxy on a nearly radial orbit. iHernquist &: Quinn 



(jl988l . Il989l ) showed that in such a case, the mass from the progenitor accumulates at the turn- 
ing points of its orbit, producing shells of high-density material at nearly constant radius on op- 
posing sides of the host galaxy. The dynamics governing the formation and shape of the shells 
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gravitating, cold, co llisionless matter forms a serie s of infinite-density peaks at successive radii. 



known as caustics. iMohayaee &: ShandarinI ()2006l ) extended this case to include warm matter 
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with various velocity dispersions: the effect of the velocity dispersion is to make the peaks of fi- 
nite width and height, so they are no longer caustics in the mathematically rigorous sense but 
retain many of the same properties, including the possibility of extremely large local density 



through self-annihilation from infall caustics in the dark matter halos of galaxies (e.g. Hoean 
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erwise determine ways in which dark matter in caustics could be detected (Duffv &; Sikivie 


2008; 



Nataraian k Sikiviel l2006l : ISikividl2003l ). Most work on the gamma-ray signal has found that caus- 
tics enhance the production from a smooth distribution by a factor of between 10 and 100. In the 
case of so-called "tidal caustics" like those seen around shell galaxies, the infall is not spherically 
symmetric, but the density is still clearly enhanced in the shells. 

Dark matter and luminous matter alike are concentrated in tidal caustics a t large distances 
frorn the bulk of the luminous matter in the host galaxy. Shell galaxies (see, e.g., iMalin k Carter 
19831 ) are beautiful examples of the extreme case of this phenomenon, where the angular momentum 
of the progenitor is nearly zero. Unfortunately, all the known shell galaxies are too far away to 
indirectly detect the dark matter in the shells: the flux of gamma rays is too attenuated, incoming 
charged particles are deflected by the Galactic magnetic field, and high-energy detectors have 
insufficient angular resolution to separate the shells from the host. However, M31 also appears to 
have shells (jFardal et al.ll2007l ). though they are not as symmetric as those in classic shell galaxies, 
and is close enough that the Fermi LAT can distinguish their position from t hat of M31's center 
(jRando et al.ll2009l ). Furthermore, an N-body model of the shells already exists (jFardal et al.ll2006l ) 
that can be used to estimate whether dark matter in them could be indirectly detected. 

N-body models of dark matter distributions have been used to estimate standard-model particle 
fluxes for indirect detection in the Galactic halo ar id the factor by wh ich dark matter substructure 
could increase the rate of pairwise annihilations (jKuhlen et al.l 120071 ). Both these quantities are 
proportional to the volume integral of the square of the dark matter density, which we will call 
the "rate" for short. The rate is estimated from an N-body representation of the dark matter 
distribution by substituting a Riemann sum for the volume integral and inferring the density 
in each Riemann volume from the N-body representation by one of many well-studied methods. 
However, neither the estimation of the square of the density rather than the density itself nor the 
choice of a suitable set of Riemann volumes has been tested. Likewise, the ability to recover the 
correct rate when the density distribution is sharply peaked has not been explored, although this 
is the scenario that would most likely lead to an observable signal and the reason that the M31 
shells are of interest. 

This paper describes tests of a number of well-known algorithms for calculating the rate and 
discusses the best algorithm to use in situations where the density gradient is large (Section [2]). We 
then present estimates, using the optimal algorithm, of both the boost factor from the M31 shells 
over the smooth distribution of dark matter in M31's halo and the rate at which gamma rays from 
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pairwise annihiliation would be seen by Fermi given likely parameters of a supersymmetric dark 
matter candidate (Section [3]). 

We find that the best way to estimate the rate from an N-body representation, whether the 
density is nearly uniform or has a large gradient, is with the simplest possible method: a nearest- 
neighbor estimator with a relatively small smoothing number to find the density, and fairly small 
constant Riemann cubes to perform the integral. This result is surprising, given that so many more 
sophisticated density estimators exist. We further find, using this result, that the largest boost 
factors from tidal debris in M31 are 2.5 percent in the most concentrated regions of the shells, and 
that the gamma rays from the debris are too few to be detected by Fermi for likely supersymmetric 
dark matter candidates: the total additional flux in gamma rays, which is model-dependent, is less 
than 7.4 x 10~^^ 7 cm~^s~^ for likely dark matter models. 

An ancillary result from our analysis of the tidal caustics, and a consequence of radial infall, 
is that the density profile of each shell and the radial spacing of the shells depend on the radial 
derivative of the gravitational potential at each shell and on the mass and size of the dwarf galaxy 
before infall. This means that information about the initial qualities of the dwarf galaxy can be 
inferred from the shells without requiring a detailed model of the potential of the host galaxy. 
Assuming that the stars in the dwarf galaxy are initially virialized (with a Maxwellian velocity 
distribution), the density profile can be fit with an analytic function whose width depends on these 
properties. We discuss further implications of this result for recently discovered shells around other 
nearby galaxies. 



2. The Optimal Estimator for High Density Contrast 



The key to the calculation presented in this paper is the estimation of the integrated squared 
density, 

r = J p'dv, (1) 

from an N-body representation of the dark matter mass distribution. The particles making up 
the N-body representation are independent observations of the mass density function p. Generally, 
proba bility density f unctions are defined as those that are everywhere positive and normalized to one 



(as in Ilzenmanll2008l . Chapter 4). The mass density function sampled by the particles in the N-body 
representation satisfies the first of these two conditions, and dividing by the total mass to get a scaled 
number density will satisfy the second. So the analysis of estimators for the probability density and 
its functionals applies equally to the problem at hand. The development and characterization of 
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Density estimators studied in the literature are divided into two classes: parametric (in which 
a particular functional form for p is assumed) and nonparametric (in which assumptions about 
the form of p are kept to a minimum). Nonparametric estimators are commonly used with data 
sets like N-body realizati ons, where the goal is usually to discover the form of p and/or calculate 
other quantities from it (jlzenmani |2008|). Among the wid e variety of nonparametric estimators 



available, nearest-neighbor estimators (IFix &: HodgesI Il95ll ) are one of the oldest and most well- 
studied varieties. The nearest-neighbor estimator uses an adaptive local smoothing length equal 
to the distance to the A^*'^ nearest particle to the location where the density is being estimated. 
The density at that point is then taken to be Ng/Vd{Ns), wh ere Vd is the volume in d dim ensions. 



cent ered on the target l ocatio n, that encloses Ng particles. iLoftsgaarden fc Quesenberrvl (jl965l ) 
and iDevrove &: Wagner] (j 19771 ) showed that nearest-neighbor density estimators converge to the 
underlying distribution at every point as the number of particles in the realization, Np, goes to 
infinity, provided that Ng/Np — )• in the same limit. They can also be considered as part of the 
larger class of adaptive kernel estimators (jMoore &: Yackellll977l') and are even more clo sely related 
when Vd is replaced by a weighted sum over the Ng particles (jMack &: Rosenblattlll979l ). However, 
because the function they return may not be normalizable, nearest-neighbor s is more suit ed to 
individual density estimates at a point than to recovery of the entire function (Ilzenmanlll99ll ). All 
the estimators we test in this work are based on either the simple nearest neighbors method or one 
using a weighted sum, although the shape of Vd varies. We describe them in detail in Appendix [Al 

The usual measure of the quality of a nearest-neighbors estimator is its root-mean-squared 
(RMS) error, 

1 



r.m.s.e. 



r 



true 



(2) 



The RMS error compares the expectation value of the estimator, in this case the rate estimator 
r, with the true value of the rate, Ftrue- For the tests in this w"ork, we used density distributions 
for which Ftrue may be calculated analytically. iBickel &: Ritovl (|l988l ) demonstrated that, given 
some constraints on the maximum slope of the underlying density distribution, the error of a one- 
dimension al integrated s quare d density estimator with a kernel of a constant size can converge 
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as Np ; iGine Sz Nickll (|2008l ) recently showed that a simple estimator of this type can be made 



adaptive using a particular rule to calculate the kernel size from the data and still converge at 
the same rate. Most of the est imators we te st in this work use adaptive kernels with a simpler 
rule than the one suggested by iGine &: Nickll for two reasons. The first is simply conceptual and 
computational simplicity: rules for choosing an optimal kernel size tend to require minimizing the 
cross-validation function (a proxy for the RMS error) of the data, which requires an optimization 
program, and the resolution convergence even with the optimal kernel chosen i n this way can still be 



— 1/2 

slower than Ap . The second is that extending t he result of IGine &: Nickll to several orthogonal 



dimensions is not trivial (jWu. Chen. &: Chenll2007l ). 



As is common in the literature, we consider the RMS error in two parts: the bias and standard 
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deviation (jLindgrenl 1 19761 ) . where 



(r.m.s.e.) 



std(r) 



(3) 



The bias, b, is the difference between the expectation value of the estimator and the true value of 
the parameter it is estimating. An unbiased estimator has 5 = 0, one for which E{T) > Ttrue has a 
positive bias, and one for which E(T) < Ttrue has a negative bias. The standard deviation indicates 
the size of the spread of individual estimates around the expectation value. For this work we scale 
the bias, standard deviation and RMS 6rror by sl fB-ctor of Ttrue; so 
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r-Tt 



EiT) 



1 



(4) 



and 



std(f) 



r 



IE 



true 



(t-E{t) 



(5) 



are consistent with Equations ([2]) and 



We used numerical experiments to assess the bias, standard deviation, and RMS error of the 
various estimators, so we must be clear about how these values are calculated numerically. For each 
experiment, 10^ random realizations of the density distribution of interest comprise one sample. 
The expectation value of a quantity is then defined as the mean of that quantity over the set of all 
random realizations. The random realizations are subject to Poisson fluctuations, so this number 
of realizations corresponds to sampling error of about one percent. We take 20 samples of the 
expectation value, so the relative error on the mean from these 20 samples is about 0.2 percent. 

The number of particles in each random realization (shown as points in an example in Figure 
[1]) is drawn from a Poisson distribution with a specified mean, denoted in the following sections as 
A'p. This precaution keeps the number of particles in the density distribution, and in the subset 
of that distribution used for the volume integral (the shaded box in Figured]), purely Poisson; the 
error associated with using a fixed number of particles depends on A'^p, so we must eliminate it if 
we wish to establish how the estimators behave as Np varies. 

The method for estimating the rate has two distinct parts: how to determine the number and 
placement of the Riemann volumes making up the sum, and how to estimate the density in each 
Riemann volume. We tested five different rate estimators that together use three different well- 
known density estimation methods and two different ways of assigning Riemann volumes (adaptive 
and constant). The rate estimators are described in detail in Appendix |A] and briefly summarized 
in Tabled! 

We first evaluate and, if possible, eliminate the bias. We constructed rate estimators using 
density estimators that have a very small or zero bias when used to estimate the density, but 
we demonstrate in this section that they do not always produce unbiased estimates of the rate 
without further correction. We wish to reduce the bias of the estimators when it is possible to 
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do so without increasing their standard deviations. Bias resulting from the statistics of Poisson 
point processes, here referred to as "Poisson bias," can be ehminated this way — analytically for 
some of our estimators, and numerically for the others — once it has been measured using random 
realizations of the uniform density distribution (Section 12. ip . We examined the RMS error of rate 
estimates for the uniform density distribution, after correcting for the Poisson bias, to separate the 
contribution of Poisson processes to the overall RMS error of each estimator from additional error 
that arises when the density is not uniform (Section 12. 2p . 

The second step in determining the best estimator is to determine the RMS error in the case 
where the density distribution has sharp features with high contrast (like shells) . Discreteness effects 
will then introduce additional bias that depends on the resolution, the smoothing number, and the 
scale of the high-contrast features (Section l2.3p . To understand when and how this bias contributes, 
we tested each estimator using random realizations of a simple caustic density distribution that has 
an analytic expression for Ftrue (Section l2.4p . Finally, we compared the RMS errors of the estimators 
for the high-contrast density distribution to determine which one to use in the calculation of the 
gamma-ray flux (Section 12. 5p . 

2.1. Eliminating Poisson Bias 

Using Poisson statistics, it is possible to construct a rate estimator r„ that is unbiased for a 
uniform density distribution; that is, one for which E{Tu) = Ptrue (Appendix[Al Equation (jASP ). It 
uses a constant Riemann volume dV for the integral and the distance to the N^^ nearest neighboring 
particle, denoted r^vs, to estimate the density. However, this estimator is not necessarily unbiased 
for non-uniform distributions. If the density is location-dependent, the integration volume dV 
must small enough to accurately sample the density gradient everywhere in the distribution. This 
choice of dV can be impractically small for distributions with high density contrast. One solution 
is to choose dV adaptively based on the local density of particles (Figure [2]) , using smaller boxes 
in higher-density regions, but this method introduces bias because the box dimensions, like rjvs, 
are then subject to Poisson statistics. We chose r„ (Appendix [Aj Equation ()A6P ) to isolate the 
contribution from choosing dV adaptively. 

Furthermore, the simple nearest-neighbor estimator is not the only option. Many other density 
estimation algorithms exist, often optimized for much better performance. However, even for these 
algorithms it is usually the case that E{fif' ^ E{n?'), especially since the bias and standard deviation 
are usually optimized for the first moment (the density) at the expense of the higher moments 
(including density-squared). We have chosen two other methods besides nearest-neighbor, both 
well-studied as density estimators, to see whether a good density estimator also makes a good rate 
estimator (Fj, Equation (|A8p . and F^, Equation ()A14p ). We also include a rate estimator based 
on one of these density estimators that has been used in the literature to calculate the rate (F^, 
Equation (lAlSP ). All the rate estimators are summarized in Tabled! 
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We calculated the expectation value and standard deviation of each rate estimator using ran- 
dom realizations of a uniform distribution, as described above. Equations ([3]) and ^ then give the 
RMS error and bias in terms of these quantities. As expected, the analytically unbiased estimator 
r^j has a numerically confirmed bias of zero (orange squares in Figure [3l). Furthermore, using an 
adaptive Riemann volume does change the bias (green diamonds in Figure [3|) . The more compli- 
cated estimators require even more substantial correction (cyan triangles, blue pentagrams, and 
purple hexagrams in Figure [3]). Interestingly, all the bias curves have the same overall shape. We 
use the multiplicative factor that transforms the biased density-squared estimator into the unbiased 
one. 



(6) 



generalized to the form 



as an ansatz for the shape of each bias curve. This parameterization evokes the introduction of an 
effective Ns, but can also adjust for the statistical effect of choosing adaptive Riemann volumes, 
and is a good fit to all the bias curves (Figure [3]) . 

Fitting E(rn)/Ttrne to a function of the form in Equation ([7]) determined 6„, the bias from 
using adaptive Riemann volumes. The results for Tf, Tg, and F^ were fit using the same function 
to determine bj, bg, and bd respectively, to identify bias from using a particular density estimation 
method. The fitted parameters for each estimator are summarized in Table El We used the fits to 
correct for the Poisson bias in all the remaining work discussed in this paper. 



2.2. Performance of Estimators on the Uniform Distribution 

All the estimators can be corrected to measure Ttmc ^or a uniform density with an accuracy of 
better than one percent, although the / and s estimators have higher-order uncorrected behavior 
(Figure H]) . The typical standard deviation of any of the estimators is larger than this residual bias 
by about an order of magnitude, so it will dominate the RMS error (compare Figure [5] to Figure H]). 
We want the estimator with the best combination of small RMS error and small Ng-. although the 
value of Ng is not as important in the uniform-density case, it limits the sensitivity of the estimator 
to small-scale density fiuctuations if the density is not uniform. We will discuss this further in 
Section 12. 3[ 

It is not clear that reducing the bias will automatically reduce the RMS error, since correcting 
for the bias can change the standard deviation of an estimator. In our case, for a given Ng, the 
bias correction simply multiplies the uncorrected estimator by a constant factor, so the expectation 
values are related by 

^(f eorr) = E{t^^,,„)/{b{Ng) + 1) (8) 
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Table 1. Various estimators of the gamma-ray emissivity T. 



Estimator 



Definition 



Long name 



See equation(s) 



<x ^^^^{hs{ri))'^dVi Epanechikov kernel density estima- IA9[ IA13| IA14I 



Unbiased nearest-neighbors estima- IA4t 
tor with uniform box size 
Unbiased nearest-neighbors estima- IA4[ IA6I 
tor with adaptive box size 
FiEstAS estimator [Ml IMl 



(jAscasibar Sz BinnevI |2005| ) with 
adaptive box size 



tor with adaptive box size 



Method from lOiemand et all (jgOOT)) IKTEl 



Note. — See Appendix El for longer descriptions of the various estimators. 



Table 2. Best-fit bias corrections. 



Estimator 


bi 


b2 


b3 


u 


1 


2 


1 


n 


1.0001 ±0.0004 


2.81 ±0.01 


-2.99 ±0.03 


f 


0.991 ± 0.004 


1.41 ±0.03 


1.41 ±0.03 


s 


0.971 ± 0.003 


-0.3 ±0.2 


4.1 ±0.1 


d 


1.00024 ±0.00005 


0.60 ± 0.01 


1.89 ±0.01 



Note. — Error ranges indicate the 95-percent confidence 
level. The form of the bias correction is given in Equation ([7]) . 
See Appendix jXl and Table [T] for descriptions of the various 
estimators. 
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and the standard deviations are related by the same factor: 



Std(f corr) = Std(f . 



uncorr 



/{b{Ns) + 1)) = std(r 



uncorr 



)/ms) + 1) 



(9) 



Prom this last expression, we see that if the uncorrected estimator underestimates Ttruej correcting 
it will increase the standard deviation. However, all the uncorrected estimators overestimate T^rue 
(Figure [3]), so correcting the bias will also reduce the standard deviation. b{Ns) + 1 is always close 
to unity and never larger than 2, so the change to the standard deviation is slight. 

For a given Ns, the RMS error of the nearest-neighbor-style estimators is generally smaller 
than that of the kernel-based estimators and decreases faster with Ng (Figure [5]) . In a kernel- 
based estimator, because each particle within the smoothing radius is individually weighted, the 
estimator must know the location of every one of the Ng particles used in the density estimate, not 
just the A'^sth one, and each weight is less than 1. Ng in the kernel-based estimators must therefore 
be much larger to get the same RMS error as in a nearest-neighbors estimator with a given Ng, 
as can be clearly seen in Figure [5l To achieve the same RMS error as r„ at Ng = 10, F^ must 
use Ng ~ 25. The RMS error of F^ converges so slowly that although it starts out with slightly 
better performance than F^, it only begins to compete with Tu at Ng > 30. Using such a large 
Ng is a built-in disadvantage for these estimators if one hopes to retain sensitivity to small-scale 
density fluctuations, and also significantly increases the computational load. For this reason we 
decided not to test the kernel-based estimators on systems with high density contrast, since the 
nearest-neighbors estimators are better suited to our needs. 

Using an adaptive Riemann volume with the spherical nearest-neighbor estimator, as is done 
in r„, increases the RMS error at low Ng. At the same time, Tf uses the same adaptive Riemann 
volume and achieves a lower RMS error. This unusual behavior may be caused by the different 
shapes of the density estimation volume (spherical) and Riemann volume (orthohedral) used in F^; 
we discuss this possibility in Section 12.41 The RMS error of F„ may be smaller relative to that of 
the other estimators in the case of high density contrast because the adaptive Riemann volumes 
can resolve the density gradient better, so we tested it, along with F^ and F/, on samples with high 
density contrast. These tests are described in the next few sections. 



In regions of high density contrast, like caustics, the maximum resolvable density is limited by 
the minimum nearest-neighbor distance expected for a given smoothing number A^^ and number of 
particles A'^p. In the limit that Ng and A'p are both large, the expectation value of the minimum 
nearest-neighbor distance scales as 



2.3. Additional Bias for Systems with High Density Contrast 





(10) 
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This scaling is derived by calculating the first order statistic of the probability distribution of 
nearest-neighbor distances for a three-dimensional Poisson point process (for a longer explanation, 
please see Appendix [B]). The scaling with A'p is valid for A'^p > 10^, but the scaling with Ng only 
approaches the asymptotic limit for values much too large to be practical (Ng > 10^). For reasonable 
values of Ng, the power-law index must be determined numerically as discussed in Appendix [Bl 

^(riv.,min) oc 7 = 0.51 ±0.06 Ap>10^ 10<A, <45 (11) 

The corresponding maximum density then scales as 

A, Ap 

nmaxOC-3 ^ ' (12) 

' Ns,niin s 

For As in the range of interest, 87 — 1 ~ 1/2. As expected, using more particles or a smaller 
smoothing number increases the sensitivity to small-scale fluctuations and the maximum density. 
The upper limit on the density introduces bias into the density estimation that also depends on 
Ap and A^ in the combination given by Equation ()12p . The estimator will perform normally as 
long as the local density is less than n^nax, but returns nmax for densities larger than nmax- Given a 
density estimator n, the undersampling-limited density estimator h^^l that incorporates this effect 
can be written 

J h E{h) < Tlmax /, ON 

The upper limit on the density changes the way the corresponding rate estimator works, since now 
the piecewise function separates regions of the Riemann sum where the density is less than 
the upper limit from regions where the density is too high to be resolved. So given a bias-free 
uniform-density estimator F, the corresponding density-limited estimator is 

f di = ti + (14) 

6 + 1 ^ ' 

where 6 is a factor with the form of Equation ([7|) and the appropriate fitted constants from Table 
[2j There is therefore an undersampling bias, in the rate estimator that depends on both Ap 
and Ng. The Ap- and A^-dependence enter two ways: in the criterion for separating the Riemann 
sum and directly in the rate calculation for one of the terms: 



Ftruc 

E{t) 



ra<Wmax _|_ ^max^>"max _ -j^ 

^true Ftrue(^ + 1) 

(15) 



Each of the first two terms is less than or equal to 1 because they are both evaluated over subsets 
of the full integration volume. Additionally, their sum must be less than or equal to 1 because 
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the limiting density is less than or equal to the density in all the Riemann volumes in that sum. 
So bus ^ 0. In the limit where the realization is fully resolved, the bias should be zero since 
Hi>nmax = 0) independent of A^^ and A'p. 

In general, the A'^p- and A'^s-dependence in Equation (llSp is complicated since some unknown 
fraction of the total volume is under-resolved. To determine the undersampling bias, we tested r„, 
r„, and Tf, including their respective corrections for Poisson bias, on N-body realizations of a one- 
dimensional caustic for which Ttrue can be calculated analytically. The caustic has an adjustable 
sharpness represented in terms of a velocity dispersion a: the smaller a is, the narrower and taller 
the peak in the density. The caustic density as a function of position and time may be expressed 
in terms of Bessel functions: 



X B 



{X - Xc) 



(16) 



with 



X>Xc ^ ' 

where Xc = l/4at is the position of the caustic and X is a modified Bessel function of the first 
kind. Table [3] briefly explains the parameters po, a, and a and the dimensions for the integration 
volume, and gives the values used in our tests where applicable. We derive this result and explain 
the parameters more fully in Appendix O Equation (I16p is given in terms of the mass density p, 
which is easily related to the number density n. When we construct random realizations of the 
caustic, we hold the normalization of the mass density po and the size of the integration volume 
V constant as A^p changes by setting the particle mass mp = poV/Np, so that realizations with 
different A^p will have the same Ftrue- 

Given the values in Table O Ttrue can be calculated by performing a single numerical integral 
(Appendix[Cl). We can use this simple model to vary the density contrast and the scale of the density 
variations, simply by generating N-body representations of the caustic described by Equation (fT6]l 
for different a. 

To make a random N-body realization of the caustic, A'^p particles are initially distributed 
uniformly in their initial three-dimensional positions q. Then a displacement function x{qx) is 
applied to the Qx coordinate to generate the sample (the set of blue points in Figure [6]). The 
general form of x{qx) is given in Equation (|C6p . The form for a particular sample is determined by 
setting the parameters a and t, which also set the location of the peak density of the caustic, Xc, 
and by setting a, the width of the normal distribution from which the random components of the 
particles' initial velocities are drawn. The values of a, t, and a we used are summarized in Table 
[3l The locations of the particles in y and z remain uniform. 
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It is important to choose the initial range of qx so that the corresponding range of x values 
defined by the mapping x{qx) is larger than the integration range in x, since otherwise the system 
will be incompletely sampled in x. In practice, we determined the range in Qx by choosing a 
range in x that is larger than the integration range and then inversely mapping it to Qx under the 
assumption that the random velocity contribution is zero (otherwise the mapping is not invertible) . 
This method does not work for random velocities comparable in magnitude to the bulk velocity, 
and we adjusted our method of calculating Ftrue to account for the incomplete sampling in these 
cases, as noted in Appendix O 

To test the estimators, we used them to calculate the rate for a set of samples with a given 
density contrast and resolution by integrating density-squared over the integration volume (shown 
as a yellow box in Figured]). The integration volume is smaller than the dimensions of the realization 
to avoid unwanted edge effects, but extends well beyond the edge of the caustic, which is the feature 
of interest. Since the undersampling bias depends on both the resolution and smoothing number, 
we varied both these parameters for each set of samples with a given a. Table H] summarizes the 
ranges and step sizes we used to explore this parameter space. Because the parameter space for 
these tests was so much larger than in the uniform-density case, we used 5000 random realizations 
at each combination of contrast and resolution, so the expected level of sampling fluctuations is 
about 1.5 percent. 



2.4. Understanding the Undersampling Bias 

We draw conclusions about the behavior of the undersampling bias with various Ng and A'p 
using the results of our tests at various levels of density contrast, represented in our model caustic 
by the parameter a. 

As is expected, using higher resolution (a larger A'^p) leads to better rate estimates of sharper 
features (Figure [7]). Even the highest-resolution realizations we tested could only resolve moderately 
sharp features. The estimators using adaptive Riemann volumes (for example, Tf shown in the 
right panel of Figure [7]) appear to require more particles to obtain the same bias as the constant- 
Riemann-volume case; this is partly due to the fact that in the adaptive scheme there is exactly 
one volume per particle, but the same number of constant Riemann volumes (about lO'*) is used 
regardless of resolution. 

We found that the undersampling bias depends strongly on Ng when the features of interest 
are marginally or under-resolved, and only weakly when they are fully resolved (Figure [HJ left 
panel). As expected, using a smaller Ng leads to a lower bias because the size of the volume used 
for density estimation scales with Ng as shown in Equation (llip , so that a larger Ng blurs out more 
small-scale features. A lower Ng also leads to a larger standard deviation of the density estimates, 
which increases the total RMS error, but this effect is small compared to the improvement in the 
bias in the under-resolved regime and only very slightly affects the fully-resolved case (Figure [HI 



- 14 - 



Table 3. Parameters for the analytic one-dimensional caustic. 



Parameter 


Value 


Dimension 


Notes 


Po 


625 


[M]/[L]3 


Mass density of the sample at t = 0. Particle mass is adjusted 
to keep the same mass density at varying resolution. 


a 


varies 


[L]/[T] 


Sharpness parameter. A smaller a makes a sharper caustic. 


a 


1/2 


1/[L][T] 


Describes the displacement function used to generate the 
caustic. See Appendix [Cl Equations (|C6P and (|C20p. 


t 


1 


[T] 


Corresponds to Xr< = 0.5. See Appendix [Cl Equation (|C25p. 




0.5 


[L] 


The rate is integrated from —Ly^z to ij/,^, where y and z are 
the dimensions parallel to the caustic. Avoids edge effects, 
as illustrated in Figure El 




[-0.5,2] 


[L] 


Limits of the rate integration in the direction perpendicular 
to the caustic, chosen so that the rate is integrated across the 
caustic. Shown in Figure [H 



Note. — These quantities are defined in more detail in Appendix O The units are given as 
dimensions only since they may be scaled as needed. 



Table 4. Parameter space for testing undersampling bias. 



Parameter 


Values Tested 


Estimator 


{u,n,f}, corrected for Poisson bias where ap- 




propriate using Equation ([7]) and the appro- 




priate coefficients in Tabled 




10 . . . 30, in steps of 2 


logio 


3 . . . 4.5, in steps of 0.25 


logio ^ 


—3 . . . 0.5, in steps of 0.25 
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right panel). 

We found that the algorithm used to determine the Riemann volumes adaptively was sensitive 
to the way in which the boundaries of the integration volume were treated. Boundaries closest to the 
sharp edge of the caustic (which in our test cases is parallel to one face of the integration volume) 



must be trimmed as described in Section 2.2 of lAscasibar &: Binneyi (120051 ) to avoid artificially 



overestimating the rate: without trimming, the Riemann volumes on the face of the caustic next to 
the boundary are artificially elongated into the region ahead of the sharp edge where the density 
is effectively zero. The density estimate in those volumes will then be artificially high because 
the density is assumed to be constant over the entire Riemann volume, leading to rate estimates 
with positive bias when the distribution is marginally resolved (Figure O solid [blue] lines) . Such 
"bleed-over" still occurs with a constant Riemann volume but is much less significant because the 
box size does not depend on the local density. 

However, applying the trimming algorithm uniformly to all the Riemann volumes with one 
or more faces on a boundary of the integration volume artificially underestimates the rate by 
significantly reducing the total integration volume on the trailing edge of the caustic, where the 
density is small but still nonzero (Figure [U dashed [cyan] lines). In this system, restricting the 
trimming only to the face nearest the caustic edge results in the correct bias behavior (Figure [U 
dot-dashed [green] lines). Because the choice of how to treat the boundaries appears to depend on 
the particular geometry of the system in question, it may be difficult to extrapolate the performance 
of estimators that use the adaptive Riemann volumes from our test case to systems with arbitrary 
geometry. In particular, it is not immediately clear how this method would extend to the shells in 
the dynamical model of M31, which have spherical edges near several boundaries of the integration 
volume. 



2.5. Performance of Estimators on the Non-uniform Distribution 

The best estimator; that is, the one with the smallest RMS error, has the best combination of 
undersampling bias near zero and small standard deviation for the smallest a. In the resolution- 
limited regime the RMS error is dominated by the bias. If the caustic is fully resolved, the bias 
is zero and the standard deviation, which is constant for a given A'p, dominates the RMS error. 
We have already established that a smaller A'^, improves performance in the under-resolved regime 
without substantially increasing the RMS error for resolved distributi ons. Figure [TOl shows that 
all three estiinators achieve the convergence rate of Np predicted by Bickel Sz Ritovl ( 1988 ) and 



Gine &: Nickll (|2008l ) for sharpnesses that are fully resolved at the highest resolution. 



Comparing the RMS errors for sharper and sharper caustics shows that all three tested esti- 
mators have very similar performance (Figure [TTI) . The nearest-neighbors estimator with constant 
Riemann volume ( [red] circles in Figure [TT]l converges slightly faster than the other two but once 
the distribution is resolved they are nearly indistinguishable from one another (Figure \TT\ inset). 
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The FiEstAS method, thanks to the space-filhng tree organizing the particles, is faster than the 
nearest-neighbors estimators, so if the distribution is known to be completely resolved (the regime 
shown in the inset of Figure [TT]) . then this method can be used without loss of performance to 
take advantage of its greater speed. However, in situations where parts of the distribution may be 
under-resolved and creating a higher-resolution realization is not possible, F^j should be used to 
take advantage of its ability to resolve slightly sharper features with fewer particles. 
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Fig. 1. — (Color online) A three-dimensional random N-body realization of the uniform distribution. 
All the points (particles in the realization) in the full volume are used to estimate the density; 
to avoid edge effects, the volume tessellation and Riemann sum are performed over the volume 
indicated by the shaded (yellow in electronic edition) box. 
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Fig. 2. — (Color online) Example of adaptive volu me tessellation referred to in the text, produced 
by the space- filling tree in the FiEstAS estimator ([Ascasibar &: Binnevll2005l ). The shading (colors 
from red to blue in electronic edition) indicates the way in which the tree divides the space (the 
spatial distribution of the inorder traversal). For clarity, not all the boxes with edges on the 
boundaries are shown. 
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Fig. 3. — (Color online as indicated in square brackets) All the rate estimators except for ([red] 
circles) require slightly different but similarly shaped bias corrections. Equation ([7]) is an adequate 
ansatz for the A'^s-dependence of the correction: an example corresponding to Equation ([6]) is shown 
as the [black] dashed line and the [green] filled points show the effect of using Equation ([6]) to correct 
Tf. Fit parameters are listed in Table [21 Table [1] explains the label abbreviations. 
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Fig. 4. — (Color online) Once estimators are corrected for Poisson bias using the fitting formula 
in Equation ([7]), the corrected versions all have a bias much smaller than the typical standard 
deviation of a few percent. The residual bias is comparable to that of r„. Table [1] explains the 
label abbreviations. 
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Fig. 5. — (Color online as indicated in square brackets) The RMS error of the corrected estimators 
is around five percent, and declines with increasing Ng. T j ([green] diamonds) and Tu ([red] circles) 
have the smallest RMS error at low Ng . Introducing an adaptive Riemann volume ( [orange] squares) 
increases the RMS error significantly at low Ng. The RMS error starts at a higher value and scales 
more slowly with Ng for the kernel-based estimators ([blue] triangles and [purple] stars) than for 
the nearest-neighbors estimators, an effect of the weighting function used to estimate the density. 
Table [1] explains the label abbreviations. 
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Fig. 6. — (Color online) An example of a random N-body realization of the analytic one-dimensional 
caustic used to test estimators for undersampling bias. The points are the locations of the particles 
in the realization. The shaded (yellow in electronic edition) box indicates the integration volume, 
positioned to include the complete caustic shape and avoid edge effects. For this sample, A'p = 10^, 
a = 0.01, and the dimensions of the box are given in Table [3l 
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Fig. 7. — (Color online) Using a larger number of particles resolves sharper caustics in all three 
cases: Tu (a), F/ (b), and r„ (not shown, but very similar to the two others). Regions colored white 
(dark or bright blue in the electronic version) are considered fully resolved because the bias is less 
than the typical standard deviation of about 2 percent. In the electronic version, different colors 
indicate the magnitude of the undersampling bias bug-, blue indicates \bug\ < 0.02, 0.05, 
yellow \bus\ ~ 0.1, and red \bus\ > 0.2. 




Fig. 8. — (Color online) Using a lower value of Ng leads to lower bias whether or not the realization 
fully resolves all the features in the underlying distribution (a, main figure and inset). The increase 
in the standard deviation, and therefore RMS error, from using a smaller Ng is small compared to 
the improvement in the bias when there are under-resolved features (b). If the realization is fully 
resolved, the increase in RMS error is detectable but extremely small (b, inset). 
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0.05 




Fig. 9. — (Color online as indicated in square brackets) The adaptive Riemann volume algorithm 
is sensitive to the treatment of the boundaries of the integration volume relative to the geometry 
of the density distribution. Thanks to the location of the caustic parallel to one boundary of the 
integration volume, the rate is overestimated (solid [blue] lines) at marginal resolutions unless the 
Riemann volumes on that boundary are trimmed to fit the face of the caustic (dot-dashed [green] 
lines), but trimming all the boundaries in the same manner underestimates the rate (dashed [cyan] 
lines) . 
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Fig. 10. — (Color online) All three estimators tested on non-uniform distributions achieve a con- 

— 1/2 

vergence rate of Np (dashed line) as predicted by prior analytical work. Tests with Ng = 10 
and logj^o*^ = —0.75 are shown. 
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Fig. 11. — (Color online as indicated in square brackets) All the estimators have very similar per- 
formance, r^j ([red] circles) converges faster in the marginally-resolved regime than the estimators 
that use adaptive Riemann volumes ([orange] squares and [green] diamonds). In the fully resolved 
regime the RMS errors of the three are nearly identical (inset). 
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3. Calculation of the Boost Factor and Gamma-ray Flux 

In this section we describe the N-body model of the M31 tidal debris (Section l3.ip and present 
derivations of formulae for the boost factor and gamma-ray flux (Section I3.2p in terms of the 
numerically estimated rate. We then describe how we used the results of the bias calibrations 
described in Section [2] to correct our numerical estimates of the boost factor and flux (Section IS.Sp . 
and present the results (Sections 13.41 and I3.6P . 



3.1. N-body Model 



For this work we use the N-body model of the tidal shell system in iFardal et al.l (|2007l ). 
The model uses a Plummer sphere as the progenitor of the tidal debris, orbiting in a static, 3- 
component representation of M31's potential: a spherical halo and bulge, and an axisymmetric 
exponential disk. To construct a model of the potential, the parameters of the halo, bulge, and 
disk were first fit to a rotation curve of M31 f rom several combine d sources, not including the 
tidal shells and the associated tidal stream (see iGeehan et al.l 120061 . for details). Then the orbit 
of the center of mass of the progenitor satellite was fit to the three-dimen sional position data and 
radial velocity measurements available for the stream (jFardal et al.ll2006l ). Finally, the mass and 
size of the proge nitor were constrai ned with an N-body model of the stream using the previously 
determined orbit (jFardal et al.ll2007l ). Dynamical friction and the response of M 31 to the merger are 
ignored since the mass ratio of the progenitor to M31 is approximately 1/500. iFardal et al.l (j2007l ) 



emphasize that the N-body model is not the result of a full exploration of this many-dimensional 
parameter space, but for this work we are only interested in the end configuration of the debris, 
which acceptably matches the available observations. 

We further make the assumption that there is a comparable mass of dark matter associated 
with the stellar tidal debris. This assumption is probably generous, since the dark matter in 
galaxies is thought to be much more diffuse than the stellar matter and much of it will have been 
stripped away by tides before the progenitor of this tidal debris even reaches the starting point of 
the simulation. However, for a first estimate we consider this assumption sufficient, though it is 
not a strict upper limit for reasons we will discuss at the end of this paper. 



3.2. Formulae 



In this section we derive expressions for the boost factor and gamma-ray flux in terms of the 
estimated rate from the N-body representation. 
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3.2.1. Boost Factor 



To calculate the boost factor, we must assume a halo model because dark matter in the shell 
can interact with dark matter in the halo. Denoting the halo dark matter with h and the shell dark 
matter with s, there are three terms in the total rate: 



/< 



K + risfdV 
= j {nl + 2nhns + nl)dV, 

= ^hh + ^sh + ^ss 

Equation (jlSp shows that the boost factor (3 depends on the halo dark matter density: 

Ttot — T/i/i 



(18) 



/3 



1 



hh 



ns(2n/j + ns)dV. 



(19) 



If Hh ^ ng, the additional emissivity from the tidal debris is dominated by the first term in equation 
(jl9p and the boost factor scales linearly with the dark matter density in both the halo and the 
shell. 

F or the halo, we use th e same density distribution that was used in the dynamical model of the 
shells (jGeehan et al.ll2006l ) with the addition of a small core with size rcore of half the size of one 
Riemann volume. The core eliminates the inf inite-density cusp at r = 0. This halo is spherically 



symmetric and of Navarro- Frenk- White form ([Navarro et al.lll996l . 119971 ): 



nh,o 



[{r + rcore) /rh][l + {r + ^core 



(20) 



where nh.o = Phfi/fnp = 3.67 x 10^ kpc^^ and r^, = 7.63 kpc are determined by iGeehan et al 
by fitting a mass model to a set of measurements of dynamical tracers of M31's halo. We divide 
the fitted mass density p^^Q by mp, the mass of the simulation particles, to get consistent number 
densities n/j and n^. 



3.2.2. Gamma-ray Flux 



We use the notation of lFornengo et al.l ()2004l . hereafter FPS) to present the calculation of the 
gamma-ray flux. The differential flux of photons in an infinitesimal band of photon energy E^^ , 
d^^/dEy, can be factored into a contribution from "particle physics" that specifies the spectrum of 
the radiation and an achromatic contribution from "cosmology" — the shape, size, density. 
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and distance of the dark matter — that sets the normahzation, as in eq. 1 of FPS: 



SUSY 



^7 "--^7 

The rate at which gamma rays would be detected by the Fermi LAT is 



IT 



d^ 



SUSY 



- — — ^cff {E^ ) dE^ 



(21) 



(22) 



where Af,Q{E^) is the effective area of the detector, E^^ is the lowest detectable energy, and 
is the dark matter mass. For the Fermi LAT, Eth is 30 MeV, and above 1 GeV the effective 
area for diffuse events is roughly independ ent of energy to w ithin about 10 percent of the mean 



value over the energy range of integration (j Rando et al. 20091). Furtherrn ore, the supersymmet 



trie 

calculations of the particle physics contribution in iGondolo et al.l (j2004l ). which we use for this 
work, take ii^th = 1 GeV. So for the remainder of this work we will use E^ii = 1 GeV and assume 
j4efj is independent of energy. With this simplification we can work with the total flux for now and 
later multiply it by ^eflf to get the detection rate in the LAT. 



Also following FPS, the particle-physics contribution is 

d^SVSY ^ ^^^^ 



dE^ 



2m^ dE^ 



(23) 



where {av) is the velocity-averaged cross section and is the yield, called in iPeirani et al 
(j2004l ) and FPS. We have moved the constant factor 1/4tt to be part of (1)'=°'^™° because it is most 
easily understood as part of the attenuation of the gamma-ray flux over distance. 



To evaluate , we use a subset of the b enchmark mode l s oflBattaglia et al.l (|2004l ) to span the 
space of supersymmetric WIMP candidates. Gondolo et al. ( 20041) ha ve calculated the gamma-ray 
yields above 1 GeV for 10 of the 12 models in lBattaglia et al.l (|2004l ): models A' through L' with 
the exception of models E' and F' . We use the values given in the table for the nur nber of photons 



i n the continuum emission times the cross section, A''^,cont. W^), listed in Table 1 of iGondolo et al 
(j2004l ). The continuum emission is not as diagnostic as the line emission at E^ = my./2, but the 
branching ratio for line emission is smaller by a factor of 10^. 



Compared to Equations (|22|) and (j23|) . we find that 

r^x c^^susY 



SUSY 



dE. 



-dE, 



7,cont 



(av) 



7 



7 



2m% 



(24) 



For the models in IGondolo et al.l (|2004l ). A''^,cont. {cyv) is in the range 10 
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-24 _3 .-1^ 



cm" 



m-y is generally given in GeV. So a useful scaling of this formula in typical units is 



SUSY 



1.54 X lO^^cm^kpc^^s^^GeV" 



7 cont 



[av) 



10-29cm-3s-i 



VlGevJ 



(25) 
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In addition to the benchmarks, we also use the most optimistic value of $SUSY fj.^^ FPS. 
According to their Figure 8, <I)SUSY < iq-8 ^^1 the models explored, and the maximum occurs 
for m-^ « 40 GeV or A^'^^cont. (cf) ~ 10~^^ cm^ s~^. We use these values to represent the most 
optimistic estimate of the flux. 

The astrophysical factor depends on the square of the dark matter mass density p^: 

= (26) 

where x,y,z indicate physical distances in a coordinate system centered on the object, and the 
integral is over the total volume of the object. The factor I/AttcP accounts for the attenuation 
in flux over the distance d from the object to the observer. The key element in calculating the 
astrophysical contribution to the flux is thus determining the integral-density-squared J p^dV . 
We calculate this quantity from the simulation results, which consist of the locations of the 
simulation particles, each with mass nip. A numerical density estimator gives the number density 
Up of the simulation particles as a function of position, which is related to the mass density py, by 
mass conservation: 

p^ = m-^riy- = rUpUp (27) 



so that 



and 



p^dV = ml nldV (28) 



2 /" 2 



= J ni^V . ^E(r.). (29) 

It is useful to rewrite the e xpression (I29D with the units and values used in the N-body repre- 
sentation. For the simulation of iFardal et al.l ( 20071 ). nip = 1.68 x 10^ Mq and d = 785 kpc. E(r) 



is calculated in units of kpc ^ . So a useful version of (j29p for this work is 

/ nir, \^ ^ d 



X 



\W*MqJ V 785 kpc 
( E{t) \ 



(30) 



The total flux of gamma rays for a given model of dark matter and density distribution, scaled 
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to typical values in the problem, is obtained by combining equations ()25p and (j30p : 



where $ 



7,0 



X 10 cm ^ s ^ 



7,0 



N^{av) 
10~^^cm^s~^ 



nir, 



VlGeV; 

d 

785 kpc 



V W^Mq 

Upc-V' 



(31) 



9.09 X 10 m ^ yr ^. The effe ctive area of the Fermi LAT 



varies between 0.7-0.85 square meter above 1 GeV (|Rando et al.ll2009l ). 



3.3. Calibration of the Rate Estimate 

To calibrate the result from the density-squared calculation, we must estimate the number 
of particles in each shell (and thereby A^p) and the a that gives the best approximation to each 
caustic's shape. Since the shells are located at the apocenters of the particles' orbits, we can 
unambiguously determine which particles are in which shell by counting the number of pericenter 
passages, Aperi, that each particle has experienced (Figure [T2|) . There are three shells and one tidal 
stream in the system, each composed of particles with a different Aperi. Examining the system 
in a projection of phase space in the r-Vr plane, shown in Figure [T3l shows the order in which 
the shells were formed. The first caustic to form corresponds to the first (outermost) winding in 
phase space with the largest apocenter. Particles making up this caustic at a given moment have 
the lowest Aperi because this caustic marks the location of the first turnaround point for bound 
particles. The most recently formed caustic has the highest number of pericenters and smallest 
apocenter distance, and has not yet been fully filled: the innermost winding is not yet complete 
because only a few particles have had time to complete three full orbits. We chose the two oldest 
shells for analysis because they contain most of the mass. We will refer to the oldest caustic, whose 
constituent particles are shown in green and have undergone two pericenter passages, as caustic 1 
or shell 1. The second caustic, whose constituent particles are shown in red and have undergone 
three pericenter passages, will be referred to as caustic 2 or shell 2. 

Since the satellite galaxy's orbit is nearly radial, the caustics are nearly spherical. However, 
slight systematic deviations from spherical symmetry, in projection along r, can slightly increase 
the perceived width of the caustic and cause a to be overestimated. In order to use our analytical 
caustic to estimate a for each shell in M31 and determine the undersampling correction to the 
result from the density estimator, we had to correct for the slight asphericity in each shell. To do 
this we determined the caustic radius of each shell in projection along the (p direction by binning 
the particles in r and (j) and finding the r-bin with the highest number of counts for each slice in (p. 
The bins were chosen as small as possible for resolution while still being able to identify the peak 
r bin for each (j). Once the peak r was obtained as a function of (f), we fit this set of points with 
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60 R 




Fig. 12. — (Color online) Sorting particles by the number of pericenter passages, A'peri (the inset 
gives the shading [color in electronic version] scheme used here) , easily identifies the main dynamical 
structures in the tidal debris. This x-y projection shows the N-body model as it would be seen 
from Earth, with x and y measured relative to M31's center and aligned with the east and north 
directions on the sky, respectively. In this projection the edge of the younger (medium gray, red 
online) shell is sharper than that of the older (dark gray, green online) shell, but both are nearly 
spherical relative to M31's center. 
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Fig. 13. — (Color online) Looking at the projected phase space of r and Vr, measured relative to 
M31's center, lets us determine the order in which the shells formed. In the self-similar model, the 
outermost caustic forms first: we see here that this corresponds to the bound particles with lowest 
A^peri (shading/colors are the same as in Figure [12]). The youngest caustic contains the particles 
that have undergone the largest number of orbits. 
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a polynomial rc(</>) and calculated dr = r — rc{(p) for each particle to correct for the asphericity of 
the shell in <j). The polynomial was of the lowest order possible, since higher orders introduce more 
spurious spread at the edges of the fitted region. In all but one case a linear fit was sufficient. The 
process was repeated with the corrected particle radii, dr, in the 9 direction to find and subtract 
drc{9). Figure [HI illustrates this process. 

Correcting for the asphericity in this manner determines the radius rc{0,(j)) of the caustic 
surface in two steps, so that 



The bin widths in r used for caustics 1 and 2 in this procedure limit the accuracy of rc{9, (/)) to 0.05 
and 0.08 kpc, respectively. However, fitting the radial profile with the analytical caustic determines 
the caustic radius more accurately, assuming that rc{9, 4>) has sufficiently corrected the asphericity. 

After correcting the caustics for asphericity, we binne d the particl es in x = r — rc{9,(l)) to 



width and fit Equation ()C37p to the density profile to get baseline parameters, then decreased the 
bin width until the fit parameters converged. The final bin width was 0.1 kpc for both caustics, 
larger than the bins used to correct for the asphericity, indicating that rc{9, cj)) adequately represents 
the caustic surface. Using smaller bins produced no appreciable change in the fitted parameters. 

Having determined the ideal bin width, we fit the model caustic of Equation ()C37p to each 
density profile. The model has a total of four parameters: the caustic width 5r, Xc, the initial phase 
space density /o, and the phase-space curvature k, that are discussed in detail in Appendix O As 
noted there, this density profile is universal for caustics created by quasi-radial infall of initially 
virialized material. We allowed all the parameters except k to vary in the fit, but checked the fitted 
values of Xc, 6r and /o by comparing them to estimates. For Xc this is trivial; it should be close to 
zero after the asphericity correction. We used mass conservation to estimate /o, since the number 
of particles in the caustic and its geometry are both known, and energy conservation to estimate 
5r, using Equation (IC36|) and about 10-20 percent of the particles in each shell just around the 
caustic peak. For reference, the covariance term in Equation (1C35P was two orders of magnitude 
smaller than the other terms. 

For each caustic, we determined k by fitting a parabola to the phase-space profile near the 
peak. The portion of data to fit was determined by the range over which the parabola was a good fit 
to the phase-space profile. As a check, we also estimated k using the local gravity at each caustic, 
usmg Equation (IC3T]) . Both are shown in the table and agree with each other. Keeping the value 
of K constant, we then fit the density profile of the caustic, comparing the fitted values of /o and 
Xc with the estimated or expected values. The results of the fitting are shown in Table [H 

As expected, the model fits the data well all the way through the peak in both cases (Figure 
[T5|) . The abrupt drop in density behind the caustic (at negative x) is an artifact of selecting the 
particles via their phase-space profile. In both cases Xc, the correction to Tc from the fit (dashed 
line), is within a few bins of the value from the asphericity correction (solid line at zero). Xc is not 




(32) 



construct a radial density profile. We used the Wand rule (jWam 



19961 ) 



to choose a starting bin 
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Caustic 1 Caustic 2 




Fig. 14. — (Color online) Each caustic surface was corrected for asphericity by fitting the position 
of the peak bin in r as a function of angle, first in the (p direction (top row) and then in the 
direction (bottom row). The points show the peak bin location and the lines indicate the best fit 
for the position of the caustic surface. In the background is the two-dimensional binned density 
map of the caustic. A linear fit was sufficient for all but the (/>-dependence of Caustic 2 (top right 
panel) which used a quadratic fit. 
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positioned exactly at the peak of the caustic because the satelhte is not totally cold (see Appendix 
O for a more detailed explanation). The widths of the two caustics are close, but not identical, 
reflecting the slight difference in the initial velocity dispersions of the particles creating them. The 
estimate of 5r for caustic 1 is much closer than that for caustic 2; this could be because the second 
caustic is at about half the distance of the first and is thus more affected by the non-spherical 
portions of the potential. The spread of energies with r in caustic 2 is certainly both larger and less 
symmetric than in caustic 1. Both estimates are slightly high thanks to the simplifying assumptions 
described in the appendix. 




x=r-rcCf (kpc) x^x-r^^e,^) (kpc) 



Fig. 15. — (Color online) The radial density profiles of caustics 1 (left upper panel; green in 
electronic edition) and 2 (right upper panel; red in electronic edition) can be fit surprisingly well 
with the functional form in Equation (116p after they have been corrected for asphericity with the 
process illustrated in Figure [T4l The insets in the upper panels show how well the function fits the 
region right around the peak of each caustic, which is the most important region for determining 
the true width. The residuals (lower panels) from the data used in the fit (shown in a darker 
shade) are evenly scattered around zero in each case, indicating that the peak's relative height is 
accurately determined. The fit parameters are given in Table [5j 

To finish calibrating the rate calculation, we compared the characteristic widths 5r and A'p 
of the two caustics with those used in the numerical experiments to determine whether the rate 
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Table 5. Density profile parameters for caustics 1 and 2. 



Parameter 



Value for caustic 1 



Value for caustic 2 



Notes 



Np 19779 6524 number of par- 

ticles in part of 
caustic used in 
fit 

</>) (kpc) 38.4 + 0.2161 + 0.140 26.0 - 1.86^ + 0.346*2 + 0.3(/> defined in 



Xc (kpc) 



0.228 ± 0.005 



0.113 + 0.009 



6r, estimated (kpc) 0.23 

6r, fitted (kpc) 0.201 +0.005 
K, estimated*^ 4.6 x 10~^ 



0.39 

0.232 +0.009 
2.4 X 10-"^ 



K, fitted'^ (4.2 + 0.2) x 10""^ (2.31 + 0.04) x 10" 

/o, estimated^ 0.37 0.096 



/o, fitted^ (0.33 + 0.01) 



(0.086 + 0.003) 



Equation §2$ 
from fit to den- 
sity profile af- 
ter asphericity 
corrections 
using Equation 
([C36|l 

from fit to den- 
sity profile 
estimated us- 
ing Equation 
(|C3T]) 

from fit to 
phase-space 
profile 
estimated 
using mass 
conservation 
from fit to den- 
sity profile 



Note. — Error ranges on fitted parameters indicate the 95 percent confidence interval. 
"^Units are kpc (km s^"'^)^^. 
^Units are kpc^'^ (km s^"*^)"^ sr"-*^. 
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estimate suffered from significant undersampling bias. The numerical experiments used t = 1 Myr 
(we can choose the units of time and length freely). Caustic 1 has 5r = 0.20 kpc and Caustic 2 has 
5r = 0.23 kpc, which correspond to test distributions with log^g a = —0.70 and log^o c = —0.63, 
respectively. According to Table O Caustic 1 has log^g A^p = 4.3 and Caustic 2 has log^g -^p = 3.8. 
From the results of the tests in Section 12.51 "we find that for logj^g Np = 4.25, the RMS error at 
logio ^ — —0.75 is 4.0 percent using Ng = 10 with the FiEstAS estimator and 3.4 percent with 
the uniform estimator. For log^^Q A^p = 3.75, the RMS error at log^^Q = —0.75 is 5.8 percent using 
A^s = 10 and the uniform estimator. All these errors are dominated by the standard deviation, not 
the undersampling bias. 



3.4. Boost factor 

Based on the fits in the previous section, we expect the shells to be fully resolved in the 
simulation. So we are free to choose any of the three estimators we tested to calculate the con- 
tribution to the total rate from interactions between shell particles, r^,,. We chose to use Tj for 
convenience, since choosing a constant Riemann volume whose size relative to the shells' thick- 
ness is consistent with our tests would require at least 10^ Riemann volumes to fill the simulation 
volume, whereas with adaptive Riemann volumes the large low-density portions require much less 
computation time. Tgh, which represents interactions between dark matter from the shell and dark 
matter in the halo, was calculated using the density estimator hj (Equation (|A7P ) to estimate n^, 
and evaluating Equation (j20p for n/^ at the same points where is estimated. The core radius was 
set to half the size of the Riemann volume enclosing the origin; in practice about 0.1 kpc. T^h was 
calculated analytically by integrating Equation (|20p over the simulation volume. We find that the 
boost factor /3 = 2.4 x 10~^ for F integrated over the entire simulation volume, independent of the 
particle physics model. 

By far the largest contribution to Thh comes from the very center of the halo. To get a 
more realistic estimate of the boost factor we recalculated it with this region excised, which would 
certainly be done for a real observation given that astrophysical gamma rays appear to come from 
the disk. We excluded a square region 1.35 kpc (0.2 degrees) on a side, centered on M31's center. 
This choice of exclusion region corresponds to about twice the resolution limit of the Fermi LAT, 
and is intended to be equivalent to excising the central 4 pixels about M31's center. This technique 
increases the boost factor by only a small amount, to /? = 0.0027 or 0.27 percent. 

Finally, we mapped the spatial variation of the boost factor by using the space-filling tree in 
the FiEstAS algorithm. As expected, the largest boost factors come from the edges of the two shells 
containing most of the mass, as shown in Figure [161 th^ maximum is 2.5 percent, independent of 
the particle physics model. 
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Fig. 16. — (Color online) The largest boost factors, shown as the lightest shades (red in the 
electronic edition), are from the edges of the two shells. The contrast in this figure is independent 
of the parameters for the particle physics model (summarized as $^u^^). 
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3.5. Astrophysical factor 



The "astrophysical factor" is the quantity ^^"^'^o defined by Equations ([26]), i^, and pO]l . 
For a given model of dark matter, comparing values of ^^^^^^o gives the relative strength of different 
mass distributions as sources of high-energy particles via self-annihilation. As described in Section 
3.41 we used to calculate for the shells, both for interactions between two dark matter 

particles in the shell material and interactions between dark matter in the shell and dark matter 
in the halo (Table [6|) . For comparison we also calculated ^^^osmo j^j, ^j^g dwarf galaxy used in the 
N-body model. The dwarf is represented as a Plummer sphere with mass Mp = 2.2 X 10^ Mq and 
scale radius b = 1.03 kpc. Integrating the density-squared of the dwarf over volume shows that 



r 



dwarf 



¥ \8mp) 



(33) 



to be used with Equation [30l We also compared these values with one calculated by IStrigari et al 



(|2007l ) from measurements of the mass and mass profile of the Ursa Minor dwarf galaxy in the 
Milky Way, scaled as if it were located in M31. Ursa Minor h as approximately the same mass 
as M31's faint satellites AndIX and AndXII (jCollins et ahlboiol ). We find that the shell signal is 
comparable in magnitude to the signal from this type of dwarf galaxy at the same distance, and two 
orders of magnitude less than the signal from the original Plummer sphere (whose mass is about 
ten times that estimated for Ursa Minor and its analogues in M31). 



Object 


^cosmo^ Gev2 kpc cm-6 


shell-shell (using Tss) 
shell-halo (using Tgh) 
Plummer dwarf (using r^warf) 
"Ursa Minor" fbased on Strieari et al. 20071 


9.9 X 10-^ 
8.5 X lO"'^ 
4.1 X 10"^ 
6.9 X lO"'^ 



Table 6: Values of the astrophysical factor 
calculated using Equation ([SUjl . 



for various configurations of the tidal debris. 



3.6. Gamma-ray Signal 



Beyond calculating F and /3, we used Equation (ISTI) t o esti mate the flux of gamma rays in the 
Fermi LAT for the various benchmarks in lOondolo et all (2004). The results are shown in Table [7] 
along with the parameters used to calculate $SUSY each case. 

By using the three-dimensional representation of the density-squared field from the shell con- 
structed with the FiEstAS estimator as a piecewise definition of the rate, we can integrate the rate 
along the line of sight and across pixels of arbitrary size such that the sum of the fiux in all the 
pixels equals the total flux in Table [71 With this technique we made two-dimensional maps of the 
expected gamma-ray emission using the most optimistic value of <I)SUSY (labeled as "Upper Limit" 
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in Table [Tj). Using this value, even the center of M31's halo is only barely detectable by Fermi, 
if the halo is shaped as we assumed for the dynamical model and the amount of dark matter in 
the dwarf galaxy is comparable to the amount of luminous matter, although the emission from 
the dark halo completely dominates over that from the shell (Figure [T71 left panel). The tidal 
structure is at least an order of magnitude too faint to be detected even if the halo component is 
fitted and removed (Figure [TTl right panel) . Because of interactions between halo and tidal dark 
matter, the signal from the tidal debris scales linearly with the mass of its progenitor as long as 
Ph > Ps (Equation (fTOj) ) so the ratio of dark matter to luminous matter in the dwarf galaxy would 
need to be several orders of magnitude larger, even after tidal stripping, for the tidal debris to be 
detectable with Fermi or for the scaling of the boost factor to become quadratic in the shell dark 
matter density. For a smaller progenitor such a ratio might be plausible, but dwarf galaxies with 
masses of 10^ Mq or higher tend to have comparable masses of luminous and d ark matter in the ir 
centers based on our understanding of the Tully-Fisher relation at those masses (|Geha et al.ll2006l ). 



o 1 ■ 

TS 0" 



-1 



-2 



1 1 r 




-2-1 1 

X (degrees) 




^10 12 

X (degrees) 



Fig. 17. — (color online) Gamma-ray emission from the dark halo, though faint, dominates over 
emission from the shell even at large radius (left panel). If <I>^^ hh is removed, the remaining emission 
from the shell is too faint to detect with Fermi even for the most optimistic parameters in the set 
of benchmarks (right panel). In these images the pixels are 0.1 degree on a side to imitate the 
approximate degree resolution of Fermi at the energy scale of interest, and the zero of the degree 
scale is centered on M31's center. 
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Table 7. Contributions to the flux of gamma rays above 1 GeV from WIMP self-annihilation, for 

various MSSM benchmarks. 



A' B' C D' G' H' r J' K' L' Upper Limit 



m^, GeV^ 


242.8 


94.9 


158.1 


212.4 


148.0 


388.4 


138.1 


309.1 


554.2 


181.0 


40 




120 


782 


195 


63.6 


1032 


86.5 


6303 


930 


7.08 X 10-* 


1.87 X 10** 


1.30 X 


^SUSYb 


3.14 


134 


12.0 


2.18 


72.7 


0.885 


510. 


15.0 


356. 


882. 


1.26 X 


^7,hh,all 


3.38 


144. 


13.0 


2.34 


78.3 


0.953 


550. 


16.2 


383. 


950. 


13530. 


^7,sh,all 


0.00780 


0.341 


0.0307 


0.00554 


0.185 


0.00225 


1.30 


0.0382 


0.906 


2.25 


32.0 


^7,ss,all 


< 10"^ 


0.00424 


< 10-3 


< 10-3 


0.0023 


< 10-3 


0.0161 


< 10-3 


0.0112 


0.0279 


0.398 


*7,addl,all 


0.00810 


0.345 


0.0310 


0.00561 


0.187 


0.00228 


1.32 


0.0387 


0.917 


2.28 


32.4 


total, all^ 


3.39 


145. 


13.0 


2.35 


78.5 


0.955 


551. 


16.2 


384. 


953. 


13560. 


^7,hh,nc 


1.92 


81.8 


7.35 


1.33 


44.4 


0.5402 


311. 


9.17 


217. 


539. 


7672. 


^7,sh,nc 


0.00509 


0.217 


0.0195 


0.00353 


0.118 


0.00144 


0.827 


0.0244 


0.577 


1.43 


20.4 


^7,ss,nc 


< 10-=* 


0.00362 


< 10-3 


< 10-3 


0.00196 


< 10-3 


0.0138 


< 10-3 


0.00961 


0.0238 


0.339 


^7,addl,nc 


0.00518 


0.221 


0.0198 


0.00359 


0.120 


0.00146 


0.841 


0.0248 


0.586 


1.46 


20.7 


^7, total, lie 


1.92 


82.0 


7.37 


1.33 


44.5 


0.542 


312. 


9.19 


218. 


540. 


7693. 



Note. — The subscripts hh, sh, and ss refer to the various terms in Equation I II8I I. N-y(av) has units 10 ^® cm3 s ^ . 'I'suSY 
has units 10-^ cm* kpc-l s-l GeV-2. All 4>.y have units lO-*"* 7 c m ^ s * . For refer ence, the Fermi point source sensitivity 
for photons with E > 100 MeV is on the order of IQ-® 7 cm^^ g-i jRando et aIJl2009^ . 

''Taken from Table 1 of lGondolo et all | |2004|) except for the rightmost column, which is based on Figure 8 of FPS as described 
in the text. 

'^Calculated analytically using Equation II25I I 

Calculated analytically using the same NFW halo as for the dynamical model. 

''Calculated by constructing a numerical estimate for risiidi in each integration volume element, then evaluating njjaio ana- 
lytically at the center of that volume element and assuming its value is constant over the entire element. 

° Calculated numerically as described in the text. 

^'S'j, addl = hs + *7, ss. 

'^•I'7, total = "I'7. hh + *7, hs + *7, ss • 

^all: integration is over entire line of sight and covers the region shown in Figure [TTl in x and y 
'nc: a central region is excluded from the calculation, as described in the text. 
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4. Conclusions 

We find that unless all the features in a given density distribution are known to be fully 
resolved, the best way to estimate the volume integral of the square of the density (the "rate") from 
an N-body realization is to use the simple nearest-neighbors estimator with a constant Riemann 
volume. If the realization completely resolves even the sharpest features, all three estimators we 
tested should agree on the result. The simplest method for estimating the rate works best for 
this problem because the other, more complicated algorithms are optimized for density estimation, 
not rate estimation, and because estimators using adaptive Riemann volumes appear to require 
slightly more particles to resolve features of a given sharpness. In any case the estimator should 
be calibrated for Poisson bias as we describe in Section 12.11 We also find that the improvement in 
the standard deviation achieved by increasing the smoothing number is smaller than the increased 
bias from blurring more small-scale structure for Ng > 10. The correct calibration of the estimator 
for Poisson bias can change the estimated result by up to 10 percent for reasonable values of Ns- 
The correct calibration for undersampling bias can change the result by a factor of 2 or more if the 
simulation is under-resolved; rather than attempt to correct for it, it is better to ensure that the 
N-body realization has sufficient resolution for the small-scale features to be resolved. 

Using a calibrated estimator and a sufficiently resolved N-body realization, we calculated the 
boost factor and signal in gamma rays from tidal debris in M31 that displays high-contrast features. 
Although we find as expected that the largest boosts come from the shell edges, they only increase 
the total signal by at most 2.5 percent over the signal from a self-consistent smooth halo. Likewise, 
the total gamma-ray flux from the shells is three orders of magnitude lower than emission from the 
dark halo, and too low to be detected by Fermi for likely dark matter parameters (Table [7]). The 
total signal is comparable to that predicted for an ultra-faint satellite of M31. 

5. Future Work 

The existence of shell features around M31 provides many avenues other than indirect detection 
for learning about the nature, dynamics, and distribution of dark matter. The very existence of the 
shells demands that the dwarf galaxy that created them must have had very low angular momentum 
relative to M31 because the pericenter distance is so small. Whereas high-angular-momentum 
systems like that of the Sagittarius dwarf galaxy in the Milky Way are useful for constraining the 
shape of dark halos because the tidal debris explores a large range in angle, low-angular-momentum 
systems like the M31 shells and giant stream probe M31's potential over a large range in radius, 
and are best suited for constraining the degeneracy between the different mass components of the 
host galaxy. They also act as a sensitive probe of the mass profile of the progenitor, since the 
combination of relatively cold initial conditions in the dwarf and a small pericenter distance acts 
as a kind of spectrometer, spreading the mass of the satellite galaxy out in space according to its 
total energy. Because the shells' relative orientations are a good limit on the projected angular 
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momentum of the progenitor, variations in the initial position and velocity of the center of mass 
are not likely to be degenerate with variations in the shape and phase space distribution of the 
debris, although this is still being tested. This makes the shapes and phase space distributions of 
the shells extremely sensitive to the initial phase space distribution of the progenitor satellite, and 
can place limits on the cuspiness of the mass profile of the dwarf. 

The analytical caustic used in this work can also be used as a model for caustics that form under 
the much more complicated equations of motion responsible for quasi-radia l gravitational infall. As 



shown in Appendix [Cl the resulting form is identical to that obtained bv iMohavaee &: Shandarin 



(j2006l ) in their analysis of those caustics with intuitive identifications of the normalization, caustic 
location and distance from the caustic surface, and requires no numerical integration to obtain 
the complete profile so it may be easily used for fitting. Although our model is less general (it 
does not predict the relative locations of caustics) it is consistent with the more general case, and 
more tractable if only the universal density profile is desired. The height and width of each caustic 
are sensitive to the initial phase space distribution of material in the caustic, while the profile 
depends on the potential of the host galaxy only through the gravitational force at the location 
each caustic — a complete mass model is not nece ssary. In light of recent discov eries of shells around 



many more nearby galaxies besides Andromeda (i Martinez-Delgado et al.ll201Cl ) , this technique may 



provide a way to constrain the properties of luminous matter in dwarf galaxies by examining the 
tidal debris they produce, as will be discussed in an upcoming paper (Sanderson, in prep.). 

Although the M31 tidal debris is probably not a candidate for indirect detection, we hope 
that our discussion of how to estimate such signals from N-body realizations will improve those 
estimates in future work. We also hope it may inspire attempts to develop optimized estimators 
for this quantity, similar to the way that optimized density estimators have been developed, since 
no optimization other than simple bias correction was applied to the estimators used in this work, 
and connect the astrophysics community with the body of statistical literature on such estima- 
tors. N-body modeling will undoubtedly prove an indispensable component of the prediction and 
interpretation of direct and indirect detections of dark matter. 
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A. Rate Estimators 

Here we describe the five rate estimators tested in this work. They are based on three types of 
density estimators (spherical and Cartesian nearest-neighbor algorithms and a spherical smoothed 
kernel method) and two methods for determining Riemann volumes (constant size and adaptive 
tessellation) . 



A.l. Nearest Neighbor 

An N-body representation of a continuous number density distribution n{x) is a Poisson point 
process with a spatially varying mean. As such, all estimators (e.g., n) of the density and its higher 
moments {"nP, n^, etc.) obey the statistics of point processes. In the case of a uniform distribution, 
these are simply the well known Poisson statistics, and the simplest estimator calculates the density 
in terms of the distance to the Nf^ particle, called the nearest neighbor distance. For a given value 
of the smoothing number Ns there is a particular nearest neighbor distance r^g for each particle, 
and the density near the particle is estimated using 

Ns 

However, if we compute E{hh) by integrating over the probability distribution describing the dis- 
tances between particles, we find that 

Eim) = (A2) 

indicating that the estimator is biased since E{hf,) ^ n. This is a result of the random fluctuations 
in the distance r^g from particle to particle, which obey Poisson statistics. Even in cases where 
the density is not uniform a similar effect is present. 

The Poisson bias of the estimator ()Aip can be easily eliminated in the case of the uniform 
distribution by noting that E{hiy) differs from n by a constant factor only. Dividing by this factor 
produces the estimator 

hu = - o , (A3) 

which has E{n) = n. 

We wish to construct a minimally biased estimator for the rate, P = J" n^dV . It is well known 
that in a Poisson distribution E{n'^), the expectation value of the square of the unbiased density 
estimator in equation ([A3]), is not equal to n^; still, an unbiased estimator for in the case of the 
uniform distribution does exist: 



Ns 



n\ = {Ns-l){Ns-2)(^-^] (A4) 
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with Ng and rjyg defined as before. Then E{n'^u) = n^- 

Using n^, we can construct an unbiased estimator for T by using a Riemann sum over Ny 
identical volumes dV to approximate the volume integral, so that 

Nv 

tu = dvY,^u,^ (A5) 
1=1 

where n'^u,i is given by evaluating Equation ()A4p at the center of subvolume i, and the total volume 
V = NydV. Using Equation dMl), E{Tu) = Ftruc- This estimator provides a useful check that the 
code is functioning properly. 

We also tested the same density estimation method with an adaptive Riemann volume, in 
which each particle occupies its own box [Ny = Np). Each particle's Riemann volume contains all 
the space closer to that particle than any other. The size of such a Riemann volume is also affected 
by Poisson statistics, so this rate estimator will not be unbiased even if Equation ()A4p is used to 
calculate the density-squared. For simplicity, we represent the Ng dependence of the additional bias 
from the adaptive box size as a prefactor, to be determined numerically, and use Equation ()A4p to 
estimate the density-squared: 



A.2. FiEstAS 

A variation on the nearest-neighbor estimator is implemented by Ascasibar and Binney (2005) 
in their algorithm FiEstAS. We refer to this estimator as /. It too uses the A^<jth nearest neighbor, 
but instead of a spherical volume considers the volume of the Cartesian box enclosing A^^particles 
when calculating the density, so that for a particle i: 

N 

nu = (A7) 

Conveniently, the construction of the tree used in calculating dVi also chooses the Riemann volume 
adaptively in the same manner as for the estimator n. Now there are two contributions to the 
bias: the Poisson bias from using (rej j)^ to estimate the density-squared and the Poisson bias from 
determining the adaptive Riemann volumes. For simplicity, we represent the A^^-dependence of 
both contributions with a single prefactor, so the rate estimator is 

1 

ff = — r^ihfifdVi (A8) 
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A. 3. Kernel-based 

We also tested two kernel-based rate estimators. Kernel-based density estimators use a weighted 
sum to smooth over the A^^ nearest particles, so that the estimated density at location x is 

n^(x) = W{xj - X, h) (A9) 
i=i 

The smoothing vector his a, generalized nearest-neighbor distance. The vector has length \xp^^ —x\. 
For a one-dimensional spherical kernel h = rjsff and the kernel function W is nonzero when \xj —x\ < 
r^. For a three-dimensional kernel, h = xn^ — x and W is nonzero when \xj — x| < \h\. 



The type of kernel chosen can have a significant effect on the bias and variance. ISharma fc Steinmetz 



(|2006l ) have tested the bias and variance of density estimators with a variety of kernels on uniform 
density distributions, and we use their notation here. The one-dimensional kernel can be written 
in the form 

W{rrh) = !^ (AlO) 
where r is the distance from the target location, u is the scaled distance 

u = r/h, (All) 



/ is the kernel normalization 



1 

7 

and Vh is the volume enclosed by the smoothing length h. 



1 

W{u)'iTTU^du, (A12) 







Sharma &: Steinmetd found that the Epanechnikov kernel 



(A13) 

I U otherwise 

has the smallest bias and variance in estimating the density in the case of a uniform distribution. 
We use this density estimator and an adaptive Riemann volume for the rate estimator s, and again 
collect the A^s-dependence of the bias in a prefactor: 

Again, 1/(1 + 6s) is the bias when using a constant Riemann volume. 



Diemand et al.l (|2007l ) used an adaptation of the kernel-based method to estimate the rate from 
simulations of the Milky Way's dark halo and halo substructure. Starting with Equation (jA14p . 
they make the substitution hs{ri)dVi = 1 (there is one particle per Riemann volume). We examine 
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this variation of the kernel-ba sed method, ref erred to as estimator d. We include a bias-correcting 



prefactor that is equal to 1 in lDiemand et al 



1 



1 + bd{Ns 



(A15) 



i=l 



From a Poisson-statistics standpoint, the substitution implicitly assumes that E{nf = E{n^), yet 
the estimator itself is linear instead of quadratic in the density. For this reason it is expected to 
behave differently than the rate estimator (|A14p . 



B. Calculation of the Minimum Nearest-neighbor Distance 

In this appendix we derive an expression for the expectation value of the minimum nearest- 
neighbor distance, i?(riv^^niin), in the case of a uniform density distribution of particles. The scaling 
of this value with the smoothing number Ng and the number A'p of particles in the simulation sub- 
sequently determines the maximum density that can be both represented in an N-body realization 
with A^p particles and calculated using the nearest-neighbor estimator with Ng nearest neighbors. 
E{fNs,min) IS the first order statistic of the estimator f^r^, which is related to the nearest-neighbors 
density estimator n by 

--(^)""''- 

We start with the PDF of the nearest-neighbor distance, 

exp(-47rnpV3) ^ 47^V ^^^--^ ( AT,np^ \ 

which can be derived from directly integrating over the joint PDF for the nearest Ng particles. The 
PDF for each particle is Poisson. To calculate the first order statistic we also need the CDF of the 
nearest neighbor distance, 

p^.Sp) = (]^73T)T - r (iVs, W/3)] = 1 - ^ ^ , (B3) 

where r(A^) and r(A^, x) are the complete and incomplete gamma functions, respectively. 

The PDF of rTv^^min is that of the first order statistic of the PDF of the nearest-neighbor 
distance: 

Substituting the expressions for the PDF and CDF of the nearest-neighbor distance. 



.._(.M/^-(^exp(-4.n.3/3)^^J \ 3 

(B5) 
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Changing variables to y = 47rn/Lt^/3 gives us a simpler expression: 



dy. (B6) 



y represents the average number of particles in a sphere of radius //. 

The expectation value of riVs,min is then determined by weighted integration over its PDF: 

using the definition of /x in terms of y and the fact that p^iVs min(y)^y ~ Pf-Ns min(^)'^A*- 

We now use two approximations for the gamma function. The first is the asymptotic expansion 
for the incomplete gamma function at large y: 

r{Ns, y) « e-y [j/^-i + O ; (B8) 

the second is Stirling's approximation for the complete gamma function at large argument N: 

V{N + 1) = iV! « N^e-^. (B9) 

Using the asymptotic expansion, we may rewrite the integral: 

^(™.„,.) ^ (ji;) f (BIO) 

Making the change of variables t = N^y, we find 



oo 



1 /3 

^(rAT^n^in) = f T^) f\ ,^ I e-*t^p(^- Wrft (Bll) 

\4TrnJ [(Ns-l)\]^pNp''^ ' Jo 
which, using the definition of the complete gamma function, is 

3 y/"^ r[iVp(iV,- l) + 4/3] / 3 [A^p(iV, -l) + l/3]! 



(B12) 



Now we use Stirling's approximation to simplify the factorials in the ratio. We define = 
Np{Ns — 1) to keep things shorter: 

{Na + l/3)\ ^ (Af„ + l/3)^"+i/3e-(^"+V3) (B13) 
[(iV, -l)!]^p \{Ns-lf^-^e-^^^-^^]^'' = {Ns-lf"e-^" (B14) 
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Substituting these two approximations back into the expression for E{f]\f^ 



we find 



1/3 



{N^ + i/3)^.+i/3e-(^.+i/3) /3(^^ _ i)x 1/3 /^^ ^ ^/3x iV.+i/3 



47rne 



In the limit ^1/3, the second term approaches e^/^, so to leading order, E( 



(B15) 



oc iV,, 



1/3. 



3(iV. - 1) 

47rn 



1/3 



(B16) 



The minimum expected distance to the particle is roughly equal to the average expected 
distance to the {Ng — 1)*'^ particle, in the limit where both Ng and A'p are much greater than 1. 

Although the criterion A'p ^ 1 is generally satisfied, we are interested in values of Ng that do 
not satisfy the criterion A^^ ^ 1: in fact, we wish to use the smallest value of possible while 
retaining a good RMS error. We must therefore estimate the scaling with A^^in our region of interest 
by tabulating values of the integral in Equation ()B7p at various values of A^ and Ap. The scaling 
in the region of interest can then be approximated by fitting the values in the region of interest to 
a power law whose index is a free parameter. We find that for 2 < log^g Ap < 5 and 10 < A^, < 50, 
the A^s-dependence roughly obeys a power law TNs,ram 

oc N] with index 7 = 0.51 ± 0.06, where 7 
depends slightly on Ap. We confirm that the scaling of E{fi\is,iam) with Apis suitably consistent 
with the prediction (Figure [T8j) . 
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Fig. 18. — The A^-dependence of E{r N^^ana) (&> thick line) is not close to the asymptotic prediction 
(dashed line) in the range of interest, and varies somewhat with Ap. We take an average slope of 
1/2 (solid thin line) for the power law index in A^. However, the Ap-dependence of -E'(rAr^^min) (b) 
is close to and approaches the asymptotic prediction (dashed line), so we use the asymptotic slope 
of —1/3 for the Ap scaling relation. 
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C. An Analytic One-dimensional Caustic 

We would like to test the bias and variance of three common density estimation schemes in 
the context of high density contrast. The three-dimensional test distribution has a one-dimensional 
caustic in (x) and is uniform in the other two dimensions (y and z). Below we present the derivation 
of Equation (I16p and discuss the calculation of F. 



C.l. Density Profile 

Following the notation of Shandarin fc: Zeldovich (jlQsd ). we define the Eulerian coordinates 
X = {x,y,z) in terms of the Lagrangian coordinates or "initial conditions" q = {qx,Qy,Qz) for a 
collisionless nongravitating ensemble: 

x{q,t) = q + Vo{q)t (CI) 

where vo is the initial velocity vector of the particles. 
Conservation of mass states that: 

p{x,t)d^x ='^pod^q (C2) 

q—^x 

where x{q, t) . The sum is taken over all values of q for which particles end up at a given x. Assuming 
a uniform starting density po, the local density p of the distribution then evolves in time as 



^(--') = Ed^ = E^^ («) 



, .\dx/dq\ A., I 

q^x q^x \^tk Qq^ 

Caustics arise when the determinant in the denominator is zero. 

To randomly sample a simple one-dimensional caustic along x in our 3D distribution, we can 
choose a uniform distribution in y and z: 

y = Qy (C4) 
z = q. (C5) 

and an extremely simple displacement rule for x that includes a constant initial velocity VQ^xiQx) with 
a small random component Vi^iQx) along the dimension of the caustic and no further interactions: 

x = qx+t {vo,x{qx) + vthiqx)) (C6) 

Vthiqx) is a random variable drawn from a distribution /(fth) for each particle in the system labeled 
by a unique qx- This "thermal velocity" is included so that F will not diverge at the caustic surface. 
Equation ()C6p is not invertible since the Vth are random. Furthermore the particle ordering at t > 
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is unknown because of these random velocities, so the density at some point may now include many 

streams, since the initial velocity determines when the particles cross the caustic. To make sure 
wc include all the streams, wc express the sum over streams as an integral over a delta function 
(dropping the subscripts on q for clarity) : 

/oo 
dq S[x-q- v{q)t] dx (C7) 
-oo 

where v{q) = vo{q) + vth{q) includes both the uniform and random parts of the initial particle 
velocity. The delta function and integral over q form a sum over all particles that arrive at the 
location x from any g at a time t. 

We choose to represent the thermal velocities with a Gaussian (Maxwell) distribution with 
one-dimensional dispersion a: 

f{vth)dvti, = ^=L=e-''th/2-^di;th (C8) 

'yth(9) is thus a function of q in the sense that for each particle, labeled by its q, a different random 
velocity is assigned. The dispersion a could easily be time-dependent, but again we take the 
simplest case and assume it is a constant. We refer to the thermal component because 

systems in thermal equilibrium have velocities described by the same distribution, so it is often 
associated with temperature since cr^ ~ T in that case. 

To make the integral in the expression for p{x) easier to evaluate, we can also replace the delta 
function with its limit definition in terms of a Gaussian: 

S[x-q- v{q)t] = lim ^^e-[^-'?-^(«)*l'/2^' (C9) 

Combining, we find 

pKh(g)](x,i) = po r dglim^i=e-[--«--(^)*lV2^^ (CIO) 
To get the distribution p{x,t) we must take the ensemble average over vth- 

POO 

p{x,t) = {p[vthiq)]ix,t)) = dvthp[vth{q)]ix,t)f(vth) (Cll) 



We reintroduce v{q) = vo{q) + vthiq) and pull out the parts of the exponent containing vth- 

/oo -| poo 1 , 

-oo V27re2 J-oo V 27rcr2 

(C12) 

For notational simplicity, we temporarily define [x — g — fo(g')t] = A. Regrouping terms, we find a 
quadratic expression in the exponent in the inner integral: 

/oo 1 2 2 f r / 1 \ tA 1 1 

!^>o 2^^"^'^'^' i_ ^^^'^ 1 - K (2^ + 2? J - l^H I 
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The integral over vt}^ can be evaluated by completing the square on the quantity in curly brackets. 
Setting 



(C14) 



we find that the inner integral becomes 

/oo 
dvth exp 
-oo 



-a vti 



(C15) 



Using the substitution u = a{vth. — tA/2e^), we can immediately evaluate the integral to be: 



I{x, t) = exp 



4e4a2 



Replacing this into the equation for the density, we find that 

1 



p{x, t) = po / dq lim 



exp 



1 



2e2 4e4a2 



Replacing with its definition and simplifying makes the limit easy to take: 



/oo 
dq lim 
-oo ^-^0 



1 



v/27r (e2 + aH^) 



-AV2(e2+o-2t2) 



and finally we are down to the last integral: 



p{x,t) 



Po 



dq e 



-[x-q-vo{q)tf/2(T^t'^ 



(C16) 



(C17) 



(C18) 



(C19) 



As (7 or t approaches zero, this expression approaches S[x — q — vo{q)t], and we recover the density 
in the case of zero dispersion. 

To perform the integral over q (that is, to sum over streams) we must choose a form for vo{q). 
The simplest form that creates a localized caustic surface is vo{q) = —ac^: 



For this form of vq we can evaluate the integral 

Iix,t)= e-[--^+°«'*]'/2-'*' 

J — oo 



(C20) 



(C21) 



in terms of Bessel functions if we make a slightly suspect change of variables. First we complete 
the square in q inside the brackets: 



I{x 



-c 



a 



dq exp<( 



1 



at AaH"^ 



+ 



1 

2^ 



(C22) 
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Then we make the substitution 



u ■ 



1 

2ai 



(C23) 



This change of variables is certainly suspect as t — t- 0, but we make it anyway, hoping that since 
we know the answer at t = we can check the result and verify that it gives po everywhere. With 
this substitution 



I{x,t) = / 

J —I 



du exp 



a I X 
2^ Wt 



1 



4a2t2 



u 



(C24) 



We now identify the caustic position Xc- In the case of perfectly cold initial conditions, the location 
of the caustic is determined by setting the denominator of Equation ()C3p to zero and solving for 
Xc, with the help of Equation (ICip . We use this same quantity to describe the "position" of the 
caustic in the warm case, since although the caustic now has a width, the peak density will still 
occur near Xc- For our choice of VQ{q) we find that this corresponds to 



Xc 

at 



1 



so that Xc = l/(4at) and 



X 

at 



1 







(x — Xc) = X 



(C25) 



(C26) 



4a2t2 

for notational simplicity when performing the integral. So now we have a factored quartic in the 
exponent: 

«^ r 9 -i 

a /- 2^2 



I{x,t) 



du exp 



2cj2 



(C27) 



Expanding this expression and pulling out a constant term leaves us with an integral that can be 
evaluated in terms of Bessel functions: 



Iix,t) 



du exp 



a 

2^2 



tU 



xa 



(C28) 



Performing the integral, replacing x with its definition, and reintroducing the prefactor, we find 
that 

'{x - Xc)^" 



p{x,t) 



Po 



a 



t3 



with 



Biu) 



X < Xc 



(C29) 



(C30) 



where X is a modified Bessel function of the first kind. Taking the limit f — )■ recovers po 
everywhere, and valida t es ou r change of variables. This form is the same as that obtained by 
Mohavaee fc: ShandarinI (j2006l ) for the density profile of a caustic with Gaussian velocity disper- 
sion, with the substitutions a^cTy — at, Ax — )■ x — Xc, and = pQ/\/ai. That is, it describes 
the same shape as a caustic formed by secondary self-similar infall from an initial population with 
Gaussian velocity dispersion, but whose normalization varies as l/\/t and whose position varies 



- 59 - 



linearly with time. In fact, this density profile is, with a few slight changes, universal for any 
initial population with a Gaussian velocity dispersion, regardless of the equation of motion of those 
particles. This is because the shape of the phase space distribution of particles near the caustic can 
always be approximated by a tilted quadratic with curvature k. For the equation of motion used 
to derive Equation ()C29p . k = at^; for motion in a gravitational potential, k depends on the radial 
gravitational force at the caustic: 

1 



2 dvl 



29{rc 



(C31) 



where g{r) is the gravitational field or radial derivative of the potential, dV/dr. 



Inspection of Equation ()C29P shows that the curvature k influences the height of the caustic 
relative to the initial density, while the product at determines its sharpness. In a more general case, 
the spreading of the particles in phase space is not linear in time and the initial velocity dispersion, 
and at can be replaced everywhere in the density profile by a generic parameter representing the 
width of the caustic, 5r. The quantity po/a = /q in the normalization of the caustic indicates that 
the physical density at time t depends on the initial phase space density, with the caveat that in cases 
where a caustic forms in one dimension but motion occurs in more than one dimension, a remains 
the one-dimensional velocity dispersion while po must be determined via mass conservation, since 
the density profile does not account for the behavior of the population in the neglected dimensions 
or for the Jacobian associated with the volume element in non-Cartesian metrics (for example, 
spherical coordinates). 

The caustic width 6r can be estimated using energy conservation. Although each particle turns 
around at apocenter where Vr = 0, the caustic itself can have a net positive velocity because its 
location at a given time depends on the relative periods of particles on neighboring orbits. Thus 
the energy at the caustic surface, assuming spherical symmetry, is: 

E = V{r) + ^v', (C32) 

For two particles on neighboring orbits, their energy difference AE is found by subtracting their 
energies, so 

AE = ^{vl-vj) + V{r2)-V{n) (C33) 

Particles at the very face of the caustic are close in r and far from r = 0, so we can write r2 = ri+Ar, 
taking Ar/ri to be much less than 1. We cannot also write V2 = vi + Av and Av/vi <C 1, since the 
average velocity of particles in the caustic may not be much larger than zero. Using the condition 
on ri, the energy difference can be rewritten in terms of the local gravity (defined earlier): 

AE=^{vj-vi)-g{ri)Ar (C34) 

In reality many particles will comprise the caustic surface; one can use this expression by taking 
its variance, realizing that r and v are correlated and assuming g{r)/Ar ^ dg/dr: 

4 = - <7(r)cov(r, v^) + g{rfal (C35) 
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The size of the covariant term is related to the amount of curvature over the region used to calculate 
the variances. It can be minimized by taking a region around the caustic small enough that the 
curvature is smaller than the radial thickness, or 5v < brj^fn. Then the covariance term can be 
ignored and the equation solved to give an estimate of the width Or'- 



Or = ^ ^-^ (C36) 

akrc) 



With these physically motivated generalizations, the density profile 



V27r V K 



(C37) 



will universally fit the projected radial density profile of any shell near its peak. Using this relation, 
we can find the location of the peak by solving dp/dr = 0. The equation reduces to a linear 
combination of various Bessel functions, which can be solved numerically to find rmax- Here we 
give the result to 10 decimal places, but more are easily available if necessary simply by increasing 
the number of iterations in the root finder. 

r-max = rc- 0.7649508674(5r (CSS) 

Substitution into the density profile gives the maximum density in terms of the physical parameters. 
Within a few percent the numerical coefficient is unity: 

/ 5r 

p„,ax = 1.0211365847/0 \ — (C39) 

V n 

The peak density occurs about one characteristic width behind the nominal caustic surface Tc 
(Figured!]). 



C.2. Rate 

Armed with an analytic expression for the density, we can now calculate the rate T for compar- 
ison with estimates from random realizations. We must integrate the density with finite limits on 
q and x, since we will use the density estimator on a finite-sized sample of the distribution. Since 
the caustic is only in the x dimension, 

T{t) = ALyL, j ^ p^{x,t)dx (C40) 

with x± the lower and upper edges of the box in the x direction, and ^Ly and the box 
dimensions in the y and z directions, respectively. Furthermore, we are now considering only 
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r 

Fig. 19. — Universal caustic form for radial infall from a population with initial Gaussian velocity 
dispersion. Formulae for the peak radius and density are given in the text. 
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particles that originated in some range [g_ , g+J , so we must also recalculate the density function 
integral, Equation ()C27p . over a finite range: 



p{x,t) 



Po 



du exp 



a 

2^ 



X + U ) 



where 



u± = q± 



1 

2^ 



(C41) 
(C42) 



Technically, we should use this expression for in the rate integral, so that 



2LyLzapl 









1 du exp 


1 X- 


.J U — 
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2^ 



[x + u 



2\2 



dx. 



with 



x± 



at 



(C43) 
(C44) 



the endpoints of the rate integral about the caustic. 



Both these integrals must be evaluated numerically, so we would like to borrow our analytic 
solution. Equation ()C29p . instead if possible. We can see that this is sufficient by examining the 
form of Equation (IC29|) . The exponential and Bessel function terms are functions of the quantity 



jx - Xc) 

2at 



Xd 



(C45) 



which compares the distance from the caustic to the characteristic length scale for dispersion, 
at. For distances much larger than at, the exponential piece approaches zero. The modified Bessel 
function /C1/4 = [l_i/4^{u) — X;^/4(n)] also has an exponential cutoff for (x — Xc) » fit, or u — )• 00, 
as seen in its asymptotic series: 



1 3 1 , . 



(C46) 



The sum of modified Bessel functions 2^±i/4 has an exponentially increasing term, shown in its 
asymptotic series: 



u S2 u-^ 
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(C47) 



that is cancelled by the exponential piece of the equation, leaving power-law convergence to zero 
for {x — Xc) ^ crt in the asymptotic series: 



e-"' (X_i/4(n)+Xi/4(n)) 



1 3 1 / 1 
n + 35^ + ^'- 



(C48) 



From examining these limits, we find that for x > Xc the density is completely negligible for 
distances more than a few times at from the caustic, and for x < Xc it becomes negligible like a 
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power law. Thus it is unsurprising that substituting the infinite-range expression for the finite-range 
one gives the right answer. 



Thus, as long as x — ^ fit at the endpoints, we can safely write 



dx (C49) 



Practically speaking, the approximation is sufficient for our needs as long as x > 500 at the 
endpoints. If this condition isn't satisfied we must do the two numerical integrals in Equation 

We change variables to x^- Inserting the form for B breaks the integral up into two parts: 

m = |£ |x.|e-24 [X_i/4 {4) +X1/4 {xi)]'dx, + |x,|e-2^^ {xj) -ly, {xl)]' dx, 

(C50) 

where we have assumed that we are integrating across the caustic so that the new endpoints, 

(C51) 

are such that Xd_ < and Xd+ > 0. This integral must be performed numerically, but it's a single 
integral rather than a two-step process. 



